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DEVELOPMENT OF AN ENGINEERING ANALYSIS OF PROGRESSIVE 
DAMAGE IN COMPOSITES DURING LOW VELOCITY IMPACT 


E. A. HUMPHREYS 


Materials Sciences Corporation 
SUMMARY 

This report describes the development and implementation of a 
methodology to predict damage initiation and growth in composite 
laminates when subjected to low velocity, low mass, lateral impact. 
The methodology incorporates a transient dynamic finite element 
analysis with composite stress and failure analyses. The procedure 
incorporates damage, as it is predicted, into the displacement solu- 
tion. Thus, the damage predicted during any time step is incorporated 
into the dynamic solution at future time steps. This coupling of 
damage and dynamic response is the heart of the computerized 
procedure. 

Utilizing a perfectly plastic impact assumption, the impact 
phenomenon is reduced to an initial velocity dynamics problem with 
the impacting mass lumped at appropriate locations. The displace- 
ment time response of the laminated plate is predicted using the 
transient analysis. Composite ply stresses and interlaminar shear 
stresses are computed based on nodal moments and forces and laminate 
models. Failure analyses are performed and appropriate elemental 
properties degraded within the displacement solution matrices. 

The analysis procedure has been utilized to simulate the impact 
of a 1.59 cm. steel sphere on an eight ply, [45/0/-45/90] g Gr/Ep 
laminate. The damage predicted included interlaminar shear fail- 
ures, laminate back face splitting and the progression of damage 
within the plane of the plate and through the thickness of the plate. 
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I . INTRODUCTION 


With the ever increasing use of laminated compos?'^ tea in struc- 
tural applications, an interesting phenomenon has becone apparent. 
This phenomenon concerns the real possibility of invisible damage 
within a composite structure caused by low velocity, low mass im- 
pacts. This type of loading environmr/nt i« most easily envisioned 
as the impact of a dropped workman's tool on a structure or, per- 
haps, runway det**is ejector^ onto a structure. 

The primary concern related to this type of loading environ- 
ment is the introduction of performance degrading damage within the 
composite structure. This damage may not be visually apparent dur- 
ing subsequent inspections of the structure and, hence, the perform- 
ance degradation will also not be immediately apparent. 

The subj<»ot of impact related phenomena has been studied by 
many investigators utilizing many different approaches. Much of 
this work has been related to ballistic type impact and, hence, is 
not applicable here. In ballistic impact, the velocities involved 
are high enough to promote large stress wave propagation effects. 

Prior work in the analysis of the low speed impact problem has 
established that it is reasonable to neglect the stress wave propa- 
gation problem and to focus on the transient structural dynamic ap- 
proach. Different approaches have appeared in the literature to 
combine contact effects with dynamic effects. The Hertz contact 
problem has been extended to the problem of dynamic contact and also 
to the problem of contact of anisotropic bodies (see ref. 1) , These 
approaches treat impact with a semi-infinite target. In the present 
case, one is concerned with a target in which the dynamic response 
of the target is important in the sense of transient structural 
motion rather than material displacement. This problem appears to 
have been addressed first by Timoshenko (ref. 2), as described by 
Goldsmith (ref. 3). Timoshenko studied the problem of the impact of 
a beam where the contact between the bodies was govern'^d by Hertz's 
law for contact deformations. Karas (ref. 4) extended the Timoshenko 


approach to the study of plate impact (see ref. 5). Moon (refs. 6 , 

7) has utilized the Hertz impact theory in combination with a Mindlin 
plate theory (ref. 8) to model a similar approach for impact of plate 
structures. This particular approach yields a nonlinear mathematical 
problem and extended numerical analysis is required to obtain solu- 
tions. Ihe procedure is sufficiently complex to motivate considera- 
tion of alternate approaches. 

Two such approaches are based on simplifications of the contact 
force analysis. In one case, it is considered that the impact takes 
place during a time period which is very short compared to the period 
of the first natural frequency. In this case, it appears reasonable 
to regard the impact as having imparted an impulse locally to the 
Plate and to then study the dynamic response of the composite plate 
carget to that impulsive loading. This approach (see ref. 9) is 
appropriate as the structural stiffness increases and the impacting 
mass deci. eases. 

Another line of approach initiated by Clebsch (ref. 10) as de- 
scribed by Goldsmith (ref. 3) assumes that upon impact, the projec- 
tile moves with the plate and that the velocity of the projectile 
becomes an initial velocity condition. Thus, the analysis is the 
structural dynamic response of the plate with the attached mass. 
McQuillen et al . (ref. 11) applied this approach to a beam. They 
minimized some of the numerical problems by considering the contact 
zone between projectilv, and target to have finite width. This ap- 
proach tended to minimize the contribution of the higher frequency 
modes and, thus, numerical procedures were more successful. How- 
ever, even with these assumptions, the work of reference 11 shows 
that modes of vibration other than the fundamental mode can be ex- 
cited by impact, particularly if the striker mass is small. This 
approach, which is expanded somewhat in references 12 and 13, is 
being utilized in the curre.>t effort. 

In the present study, the primary emphasis has focused on the 
prediction of d.amage initiation and propagation during the impact 
event and subsequent dynamic plate response. The difficulty posed 
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here ie that any induced damage will alter the plate etifinese lo- 
cally and, hence, affect the aubrequent dyneunic response. Thus, 
the closed form analytical approach taken in references 11-13 is 
not applicable. 

The approach taken has included the use of a transltint dynamic 
finite element code, modified for composites analysis. The code 
selected for this purpose was SAP IV (ref. 14). The modifications 
made within the finite element code have allowed for the computation 
of composite laminate properties, prediction of layer an<^ interlami- 
nar stresses, failure analysis and incorporation of predicted damage 
into the subsequent dynamic response. 

The finite element method is well suited for this type of analy- 
sis since spacial variations in material properties are easily incor- 
porated. Thus, local damage can be included without affecting the 
stiffness of adjacent material. 

The computerized analysis procedure, CLIP (Composite Laminate 
Impact Program), has been utilized to predict the dynamic response of 
a laminated plate subjected to low velocity impact. The predictions 
included damage initiation and growth during time of the dynamic re- 
sponse. A description of the code and users guide are included as 
an appendix to this report. 
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II, ANALYTICAL METHODOLQGy 

The analysis of laminated composites subjected to impact load- 
ings consists of several distinct procedures. The dynamic response 
of the system is required. Coupled with this are laminate stress 
and failure analyses. Additionally, predic;.ed damage must be incor- 
porated back into the dynamic response predictions. The analysis 
operates in a fashion where each procedure is dependent on the others. 
The dynamic response affects the stresses. The stresses affect the 
failure modes and locations. The failure modes and locations in 
turn affect the dynamic response. 

Each of these are is and their effects upon the other procedures 
are discussed here. 

DYNAMIC ANALYSIS 

The dynamic analysis routines used in developing the CLIP code 
were taken from SAP IV. The routines which were utilized include an 
anisotropic thin shell finite element, capable of bending and mem- 
brane loading, and a time integration scheme which integrates the 
equations of motion to predict the time dependent structural response. 

Thin Shell Finite Element 


The finite element used in the analysis combines a bending 
element with a linear curvature field and a membrane element with 
a constant strain field. The two elements are combined in such a 
fashion that bending-extensional coupling cannot be modeled. Thus, 
it is required that laminates modeled be mid-plane symmetric. 

The element as formulated in SAP IV used the same material prop- 
erties for bending and extension. Thi“; has been modified such that 
nonhomogeneous materials can be properly modeled with different prop- 
erties in extension and bending. The shell element does have the 
capability to model both extensional-shear coupling and bending- 
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twisting coupling as it utilizas a full plans stress stiffness ma- 
trix in its formulation and, as such, iz well suited for composites 
analysis. 

For dynamic analysis, it is required that a mass matrix also be 
formulated. The method used in formulating the thin shell element 
consists of generating a lumped mass vector. Hence, only diagonal 
mass terms are formed. Additionally, no rotary inertial terms are 
included. Thus, the mass is distributed' at translationary degrees 
of freedom only. 

Since dynamic responses typically may include the effects of 
material damping, some provision must be made for these effects in 
an analysis. The method used in SAP IV consists of Raleigh damping. 
This is a convenient formulation as it does not require the formula- 
tion of elemental damping matrices. The effects of damping are in- 
cluded at the global mass and stiffness level and, hence, will be 
discussed in the next- section. 

Time Integration Routines 

The transient time analysis is performed using the Wilson-0 
method (refs. 14, 15). This procedure is a modification of the lin- 
ear acceleration scheme. The method is unconditionally numerically 
stable for any choice of time step. The results of the analysis are 
dependent on the time step, however, due to numerical damping. 

The effects of numerical damping can be overcome by judicious 
time step selection. The approach involves determining the highest 
structural frequency of interest and then selecting a time step which 
is a small fraction of the period of this frequency. In reference 
15, it is shown that for the Wilson-0 method an amplitude decay of 
one percent per cycle can be expected for a time step to bending 
mode period ratio of approximately 0.045. 

The scheme used for material damping in the analysis consists 
of Raleigh damping. In this method, the damping is assumed to be 
proportional to both the stiffness and mass matrices. This is con- 
venient since no elemental damping matrices are required. Also, 
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since typically the structure to be analysed with CLIP will be of 
one material, this type of damping is most suitable. 

Initial Conditions 


The impact analysis performed by CLIP assumes a perfectly plas- 
tic impact. This implies that the impacting mass attaches to the 
plate and remains in contact. This type of impact also implies a 
conservation of momentum within the system. 

These effects are incorporated by lumping the impact mass at 
specified nodes on the plate structure. An initial velocity is then 
applied at these nodes. The velocity used is scaled such that the 
product of impact mass and impactor velocity is equal to the product 
of impact mass plus nodal masses and the applied initial velocity, 
thereby conserving the impact momentum. 

There are provisions which have been added to allow for the im- 
pacting mass to rebound from the plate structure. The contact force 
between the impactor and plate is computed based on the product of 
the accelerations of the impacted nodes and the impactor mass . When 
this force becomes tensile, the mass detaches and the analysis be- 
comes a free vibration problem. 

COMPOSITE AND STRESS ANALYSIS 

The purpose of the composites analysis is two-fold. First, lami- 
nate properties must be generated for input to the thin shell finite 
element and second, composite inter- and intralaminar stresses must 
be predicted for use with a failure analysis. The generation of 
properties and in-plane stress analysis of undamaged finite elements 
directly follows classical lamination theory and as such will not he 
detailed here. The areas which need explanation include the genera- 
tion of stress, moment, and transverse shear resultants from the nodal 
forces and moments, and the prediction of transverse (interlaminar) 
shear stresses from the shear resultants. 


6 



The etreee recovery portione of the thin fhell finite elei&ent 
in SAP IV are designed to generate both stresses and moment resul- 
tants within a triangular or arbitrary quadrilateral element. The 
capability for predicting t;ransverse shear resultants did not 
exist, how&ver. In order to include this capability and minimize 
redundancy within the CLIP code, the SAP IV st'^ess recovery routines 
were removed. The shell element procedures were then nK>dlfied to 
produce nodal forces and moments directly from the elemental stiff- 
ness matrices. By using equilibrium considerations, stress, moment 
and transverse shear resultants are generated directly from these 
nodal forces and moments. These considerations are detailed in Ap- 
pendix A. 

As a consequence of the required removal of the SAP IV stress 
recovery routines, and the inclusion of the alternate method used, 
all shell elements within the CLIP code are required to be rectangu- 
lar. Additionally, all elements must be aligned with the global 
coordinate system. 

It is still possible to use nonrectangular quadrilateral ele- 
ments and obtain the correct displacement response but the subsequent 
stress analysis would be in error. 

The prediction of transverse (interlaminar) shear stress follows 
closely the method used in classical strength of materials texts for 
beam bending shear stresses. The finite element models used do not 
include shear deformation and, therefore, the transverse shear re- 
sultants are determined through equilibrium. This method produces 
satisfactory results if the thickness of the plate is small and the 
bending response dominates the overall deflections. 

The method used for interlaminar shear stress calculation in 
both intact and damaged elements is detailed in Appendix B. Pre- 
dicting interlaminar shear stresses in damaged elements is consider- 
a!:ly more difficult and will be discussed in a later section. 

The calculation of interlaminar normal stresses within the lami- 
nated plate is possible through the use of stress equilibrium. How- 
ever, these stresses should have significant magnitude only directly 
under the impacting mass. Thus, the effort required to do a rigorous 
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analysis was deemed too extreme. As a first approximation, inter- 
laminar normal stresses are computed directly under the impact site 
as a linear function of the thickness. The contact force is dis- 
tributed over the affected elements and scaled such that the stress 
on the back face is zero. 

Failure Analysis 

As the inter- and intralaminar composite stresses are computed, 
it is necessary to evaluate whether or not failure has occurred. 

This evaluation must predict both the presence of failure and the 
type or mode of failure. The mode of failure is required such that 
the specific elastic properties affected can be modified without af- 
fecting the other material parameters. 

The failure criteria selected for use in the CLIP code are listed 
in table 1. These criteria are applied ply by ply for in-plane fail- 
ure analysis, and interface by interface for interlaminar failure 
analysis. The components of stress used in the failure analysis are 
described in figure 1. 

Incorporation of Predicted Damage 

The primary goal of the current study is to track damage accumu- 
lation and growth throughout the impact event. This involves calcu- 
lating damage at selected time steps and incorporating the effects of 
this damage back into the dynamic solution. The failure criteria de- 
scribed in the previous section are used to calculate the modes of 
failure within the shell finite elements. When the modes of damage 
are found in an element, its stiffness matrix is reformed and the 
new element is incorporated into the dynamic analysis. This proce- 
dure is accomplished by subtracting out the original element from 
the unfactored global stiffness matrix and adding the new, damaged 
element. This, of course, requires that the global stiffness matrix 
be decomposed again. 




Ply Damage 


When the damage mode predicted consists of ply damage, it is 
necessary to compute new laminate stiffnesses for bending, extension 
and bending-extensional coupling. The new bending and extensional 
properties are input back into the shell finite element routines for 
formulating the new elemental stiffness matrix. The bending- 
extensional coupling terms are used only in the stress recovery in- 
formation since, as was mentioned before, this type of material be- 
havior cannot be modeled in the displacement solution. 

The way in which the laminate properties are changed depends 
upon the type of ply damage. If the ply damage includes fiber fail- 
ure, then the entire ply is removed from the laminate model. If, 
however, matrix failure is the only failure mode, it is assumed that 
the fiber can still carry load. Thus, only the transverse properties 
of the ply are removed. These laminate modifications apply to the 
specific element under consideration. Thus, if only one element 
sustains damage, the other elements are not. affected. 

When an interior ply fails under the dynamic loading, it is 
assumed that the constraint of the adjacent material will cause the 
strain distribution to remain linear through the laminate thickness. 
Because of this, it is not necessary to compute properties for two 
or more separate laminates for the damaged shell elements. This 
type of damage does, however, pose some difficulties in finding the 
interlaminar stresses. The methodology used for predicting inter- 
laminar shear stresses is described in Appendix B. 

Within the CLIP code, the types and locations of damage are 
saved for each damaged element. This information is then used to 
evaluate subsequent damage modes. Under certain circumstances, it 
might be possible for the numerical procedure to predict the same 
damage mode and location more than once. This is physically unreal- 
istic though, and if the CLIP code senses this, the analysis is 
stopped . 


9 



Interlaminar Damage 


When interlaminar damage is predicted due to the dynamic load- 
ing, the procedure used is different from that for ply damage. Be- 
cause of the assumption that adjacent material restrains the curv^ 
ature and strain distribution when local damage is present, the only 
effect of a delamination is an adjustment to the transverse shear 
stress field. The shell finite element does not consider shear 
deformation and, therefore, the increased shear deformation asso- 
ciated with a crack tip singularity cannot be modeled. The 
introduction of interlaminar delaminations does not require a re- 
formulation of the elemental or global stiffness matrices. A full 
description of the assumptions implicit in this method, the effects 
on the shear stress distribution, and the formulation of the shear 
stress calculations is contained in Appendix B. 

SUMMARY OF ANALYSIS METHODOLOGY 

The analysis methodology described previously including Appen- 
dices A and B detail the approach taken for analyzing low velocity 
impact of composite plates. The assumptions implicit in this metho- 
dology, as well as the capabilities and limitations of the analysis, 
are summarized here for clarity. 

Analysis Assumptions 
1. Small Deflection Analysis 

The displacements of the plate structure are sufficiently 
small such that the original geometry is applicable throughout 
the analysis. 
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2. Displacement Fields 


Only bending and membrane displacements are modeled. 

The effects of shear deformation are insignificant with re- 
spect to bending deformations. Typically, this condition is 
satisfied if the plate thickness is small compared to the other 
plate dimensions. 

3. Material Properties 

Static material constants and strengths are i.sed. The 
effects of time-dependent material properties are insignifi- 
cant in a realistic structural laminate. This is primarily 
due to the presence of fibers in many directions. 

4. Superposition of Static and Dynamic Displacements 

Static pre-stress displacement fields are added directly 
to the dynamic displacements. No stability analysis is per- 
formed as a consequence of (1) . 

5. Laminate Modeling 

Bending-extensional coupling is not included in the 
analysis. This effectively requires that all laminates 
modeled be mid-plane symmetric. The assumption is made that 
damage induced bending-extensional coupling will be highly 
localized and. therefore, negligible. 

6. Perfectly Plastic Impact 

The impact involves a complete momentum transfer from 
the impacting mass to the composite plate structure. This 
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also implies that contact effects (Hertzian Contact) and 
wave propagation effects are negligible. Additionally, 
since contact effects are not modeled, the impact mass 
can only rebound due to plate dynamics. 

7. Interlaminar Normal Stresses 

Interlaminar normal stresses are computed only di- 
rectly under the impacting mass. A linear approximation 
is used through the thickness of the plate such that the 
back surface has zero stress and the impacted surface 
carries the contact force distributed over all affected 
elements. 


Damage Modes 


8. Matrix Ply Damage 


Damage produced by combined ^ 22 ' ^12 fields. 

Matrix dominated ply moduli E 2 and 

in the affected ply within the affected finite element. 


9. Fiber Ply Damage 

Damage produced by stresses. All lamina moduli 
are set to zero in the affected ply within the affected 
finite element. 


10. Interlaminar Delamination 


No stiffness effects as 
Increased shear deformations 
shear stress distribution at 


a direct consequence of (2) . 
associated with a singular 
the crack tip are not modeled. 
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Damage Stress Effecta 


11. Matrix Ply Damage 

Lamina stresses O 22 0^2 automatically set to zero 

as a consequence of (8). Additionally, interlaminar stresses 
°Kz *^yz zero within the damaged ply. 

12. Fiber Ply Damage 


Lamina stresses O 22 end 0j^2 automatically set to 

zero as a consequence of (9) . Additionally, interlaminar 
stresses o„_ and o„_ are set to zero within the damaged ply. 

X 2 y Z 


13. Interlaminar Delamination 


Interlaminar stress components and are set to 
zero at the affected ply interface within the affected finite 
element. 


Damage Propagation 


14. Ply Damage 

Damage predicted in plies propagated due to load redis- 
tribution. The load distribution is changed because of local 
stiffness changes (8) , (9) . 

15. Interlaminar Delaraination 

Delaminations do not propagate to adjacent elements. De- 
laminations may occur in adjacent elements due to shear force 
distributions but these distributions are not changed as a 
result of interlaminar damage (10) . 
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III. LOW VELOCITY IMPACT ANALYSIS 


The analysis methodology described in the previous sections has 
been utilized to predict the response of a clamped rectangular lami- 
nated composite plate subjected to a low velocity, low mass impact. 
The analyses performed have provided information relating to the 
modes and locations of impact induced damage, the displacement re- 
sponse of the composite plate and the contact force between the im- 
pacting mass and the plate structure. This information is generated 
at various times throughout the duration of the impact event. 

In the course of predicting the response of the plate, it was 
required to make several computer analyses. These included solutions 
with the entire impact mass lumped at a single, central node and so- 
lutions with the impacting mass distributed over a group of centrally 
located nodes. 

The laminate configuration, impact parameters, finite element 
models and the results of the various analyses are described and 
discussed here. 

ANALYSIS PARAMETERS 

Before discussing the various analyses performed, the various 
parameters relating to the materials, finite element models and im- 
pact conditions need to be described. The analyses performed were 
part of an effort to model one of many impact experiments which have 
been performed at NASA Langley. Hence, the materials and configura- 
tions correspond to this experiment. 

Laminate Configuration 

The laminate selected for the impact analysis was an eight ply, 
quasi-isotropic configuration. The stacking sequence analysed was 
[45/0/-45/901 s . The material properties used correspond to a T300/ 
5208 Gr/Ep system and are listed in table 2. The unidirectional 
properties used are typical static data and were taken from refer- 
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ence 16. It should be noted that variation in these properties 
could greatly affect the solution. Therefore, a design application 
of the analysis would require the characterization of the material 
under consideration to avoid the wide variation in material data re- 
ported in the literature. 

The ply thickness used was 0 01321 cm. yielding a total laninate 
thickness of 0.1056 cm. This laminate configuration exhibits bend- 
ing-twisting coupling which must be taken into account in the finite 
element model. 


Finite Element Models 


The finite element model used for most of the analyses made is 
shown in figure 2. The model encompasses the entire plate structure 
as required by the bending -twisting coupling present in the laminate 
to be analyzed. The model dimensions are 10.16 cm. by 15.24 cm. 

In figure 2, the shaded area represents the elements which were 
selected for stress analysis. The CLIP code has been developed such 
that only specified elements have stress calculations performed (see 
Appendix C) . This was done in an effort to maximize computational 
economy . 

The central section of the plate ivas chosen for stress calcula- 
tions as this is the area where any impact induced damage should oc- 
cur. More elements could have been selected for stress analysis but 
this was deemed unnecessary since the impact site was directly at 
the center of the plate. 

The model was developed to represent a clamped plate and as 
such, the edge nodes are constrained against all rotations and dis- 
placements. The model contains 704 elements, 759 nodes and 1953 
active degrees of freedom. The number of active degrees of freedom 
is minimized by constraining all in-plane displacements. 

In the finite element model, the maximum element aspect ratio 
is eight. Within the area where stress calculations are performed, 
the maximum element aspect ratio is 4 . 
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In order to verify the validity of this model, a static solu- 
tion was run. The loading consisted of a uniform pressure distrib- 
uted over the entire plate surface. The results of this solution 
were then compared with a one-term Ritz solution taken from refer- 
ence 17. For this comparison, it was necessary to assume that the 
quasi-isotropic t45/0/-45/90 J g laminate behaves as a specially or- 
thotropic laminate. Hence, the bending- twisting coupling terms, 

Di 6» D26 are ignored. Even with this approximation, the two solu- 
tions compare within 3.3%. Hence, the model has been shown to be 
valid. 

During the course of performing the impact analyses, it was 
necessary to develop another finite element model. This model is 
shown in figure 3. The only difference between the two models is 
the removal of four elements and one node at the center of the plate 
in the second model . The rationale for developing this second model 
will be discussed later in this report. 

Impact Parameters 

The impact mass and velocity modeled in the analyses are listed 
in table 3. All of the analyses used these parameters. The vari- 
ous analyses did in some cases represent different distributions of 
the impact mass on the plate, however. The different distributions 
used are depicted in figure 4. Each of the sections shown are repre- 
sentative of the geometric center of the plate. Obviously, the 
third impact mass distribution shown in figure 4 was utilized with 
the second finite element mesh which had the four central elements 
removed . 


Integration Time Step 

The selection of the integration time step is one of the more 
critical steps in defining the impact model. The time step selected 
must be small enough to adequately determine the response of all 
critical bending modes while not requiring an excessive number of 
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time Bteps to investigate the impact event. The approach is to com- 
pute the periods of the plate natural frequencies* and determine which 
are critical. 

An analysis was made using a solution given in reference 17 to 
determin^^ the natural frequencies for the clamped plate in question. 

As before, it was necessary to model the plate without the bending- 
twisting coupling terms. This implies that all frequencies calcula- 
ted are higher than actual since neglecting the coupling terms ef- 
fectively increases the plate bending stiffnesses. Additionally, 
the predictions were made without including the effects of impact 
induced damage. Thus, the stiffnesses were again higher in the 
predictions than can be expected in the impact analyses. Hence, the 
computed frequencies can be expected to be somewhat higher than the 
actual natural frequencies 

Based on these calculations, a time step of 1 ysec was selected 
for the initial impact analysis. This corresponds to approximately 
1/19 of the 10-10 bending mode period. The results of the initial 
impact analysis indicate thax. the 1 ysec time step was considerably 
smaller than required. The bending shapes of the plate did not in- 
volve frequencies this large and, hence, for the remaining analyses, 
a time step of 2.5 ysec was used with no apparent degradation of the 
results. 

In conjunction with the time step selection, it must be deter- 
mined how often to perform stress calculations. In each of the 
analyses made where stress calculations were included, the frequency 
of stress calculation was every five time steps. In terms of com- 
puting damage growth, the optimum frequency of stress calculation 
would be to compute them every time step. Considerations of the 
cost and time involved precluded this, however. 

IMPACT ANALYSES 

In order to predict the response of the laminated plate described 
when subjected to the impact conditions also described, it was neces- 
sary to perform four separate analyses. The first three analyses 
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were performed with the finite element model shown in figure 2. 

The last analysis was !?jide using the model in figure 3. 

The first analysis was made without the inclusion of stress 
calculations. This analysis was performed in order to verify 
the time step selection and determine the duration of the impact 
event. (See Appendix C.) 

The second and third analyses were full analyses including 
stress analysis. In the second analysis, the impact mass was lumped 
at the center node of the finite element model. The results of this 
analysis prompted the distributed impact mass utilized in the third 
analysis. In both of these solutions, the damage computed was so 
extensive that computerized procedure terminated the solution. 

The fourth solution was made in an attempt to model the dam- 
age growth beyond the point at which the computer program had 
terminated in the second and third solutions. Hence, the removal 
of the four central elements in figure 3. 

Convenient groupings of the damage predictions of these analy- 
ses can be found in table 4. The information in table 4 will aid 
in comparisons of the growth of damage between solutions and in 
comparisons of the types of damage within each solution. 

Displacement Solution 

The first analysis performed was simply a dynamic displacement 
response solution. The impact mass was lumped at the center plate 
node. The integration time step was 1.0 usee. 

The displacement response of the laminated plate plotted 
through the center of the plate, along the axis corresponding to 
the smaller dimension of the plate, is shown in figures 5 and 6. 

The displacements in figure S represent the very short time response 
of the plate while figure 6 depicts the longer time response. 

Comparing the two figures demonstrates a fundamental difference 
between the early displacement fields (fig. 5) and later displace- 
ment fields (fig. 6) . At very early times in the impact event, 
the displacements can be characterized as a local phenomenon. At 
T =2.5 X 10*^ seconds, the major displacement response can be seen 
to exist at the center of the plate, with the rest of the structure 
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remaining nearly motionless. As time progresses, the extent of the 
plate with significant displacements can be seen to be progressing 
outwards to the edges of the plate. 

The progress of outward spreading of the major displacement 
response is complete at T 2.0 x 10~^ seconds, as can be seen in 
figure All of the displacement fields in figure 6 can be char- 
acterized as a predominance of the third bending mode shape. The 
third mode rather than the first mode is excited due to the clamped 
boundary conditions. 

Another aspect of the mode shapes excited relates to the inte- 
gration time step selected. Previously it was stated that the 
1 Msec, time step was sufficient to characterize the 10-10 bending 
mode shape of the impacted plate. It is quite obvious from figure 
5 that the mode shapes present do not approach the 10-10 mode. 

Hence, for all remaining analyses, the time step was increased to 
2.5 Msec. 

In setting up the analysis, two unfortunate situations were in- 
cluded. First, the printing of diJ^piacements was set up in such a 
fashion that it was not possible to observe the displacement fields 
along the longer dimension of the plate. This is apparently of lit- 
tle consequence since the displacements of figures 5 and 6 have pro- 
vided sufficient information. The second problem involves the 
amount of time required for the duration of the Amppct event. The 
solution was set up with five hundred time steps. Hence, the total 
time allowed was 5 x 10“4 seconds. At the end of this time, and 
therefore the end of the solution, the impact maas had not rebounded 
and the maximum displacements had not yet been reached. This caused 
some difficulty since part of the rationale for performing this solu 
tion involved determining the time duration of the impact event. 

In order to make an estimate of the duration of the impact, it 
was necessary to consider the contact force between the impact mass 
and the plate. In figure 7, a plot of the force-time response of 
the impact event is shown. The force appears to be somewhat erratic 
as it is computed as the product of ;he impact mass and the accelera 
tion of the plate node where the mass is lumped. The accelerations 
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of the node in queetion reflect considerably higher frequencies than 
the displacements and, hence, the contact force calculated also con- 
tains the high bending mode frequencies. 

EvalUi'itions of the minimum val'te of the forces shown in figure 

7 indicate tha*- the impact mass was close to rebr nding when the 

analysis terminated. Additionally, a displacement analysis made 

with a less refined finite element model indicated that the maximum 

displacement value for these Impact conditions was nearly attained 

at 0.5 msec. Therefore, for the remainder of the impact analyses, 

the maximum time allowed for the solution was increased to 8.75 x 
-4 

10 seconds. This corresponds to 350 time steps at the new inte- 
gration time step of 2.5 psec. 

Stress Solution, Mass Lumped at One Node 

Using the new time step described, a full impact analysis with 
stress solution was performed. This solution progressed until two 
elemental laminate stiffness matrices became singular due to damage 
in all plies at the end of thirty time steps. 

The displacement time response through the center of the plate 
along the shorter dimension of the plate is shown in figure 8. Com- 
paring these displacements with those shown in figure 5 demonstrates 
little, if any, difference between the two solutions. Since in 
figure 8 considerable damage is present, especially at 7.5 x 10“^ 
seconds, one would expect significant .Ilf ferences to be present. 

The reason that these differences are not present is simply that 
the response shown in both figures 5 and 8 is primarily inertial. 

Not enough time has passed for the effects of the induced damage to 
significantly affect the solution. Had the stress solutions contin- 
ued for a longer time, the differences would have become significant. 

In setting up this solution, sufficient displacement printing 
was specified such that the response of the plate along the longer 
axis could be observed. The displacement time response through the 
center of the plate along the longer dimension of the plate is shown 
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in figure 9. The displacement fields along the long axis of the 
plate show considerable similarity to those along the shorter axis. 
The primary difference is that the changes in slope of the curves 
are slightly more gradual along the longer axis. This was expected 
since the plate is longer and, hence, more flexible in this direction. 

In both figures 8 and 9, the most striking result is the limited 
area of the plate with significant displacement response during the 
early moments of the impact event. As time progresses, the displace^ 
ment "waves" move outward until the entire plate is affected and one 
would expect high moments and shears corresponding to the compressed 
displacement fields. 

The moment resultants through the plate center along the shorter 
and longer plate axes are shown in figures 10 and 11 respectively. 

The moment resultants correspond to T ■ 5.0 x 10"^ seconds and, 
hence, the middle displacement fields in figures 8 and 9. The moment 
resultants depicted in figures 10 and 11 represent averages of the 
elements on either side of the plate centerlines. The specific 
values vary slightly on either side of the centerlines due to the 
bending/twisting coupling described previously. This coupling also 
produces twisting moment resultants. The twisting moments are not 
shown, however, since they are of such small magnitude. The maxi- 
mum twisting moment predicted was less than 8% of the peak value. 

The moment resultants shown in figures 10 and 11 demonstrate 
that the major effect of the impact is highly localized at 5,0 x 10“^ 
seconds. This is in agreement with the displacement fields depicted 
in figures 8 and 9 . 

The transverse shear resultants corresponding to the moment 
resultants are shown in figures 12 and 13, Once again, the major 
loading occurs in a local region surrounding the impact site. This 
is due, of course, to the large moment resultant gradients at the 
plate center (figs. 10 and 11). 

The transverse shear resultants are plotted in the same fashion 
as the moment resultants. The figures represent averages of the 
two elements on either side of the plate centerlines. The shear 
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resultants in individual elements also contain the effects of the 
bending^twisting coupling and, therefore, are slightly different on 
either side of the plate centerlines. 

The moment and transverse shear resultants depicted in figures 
10>13 are sufficiently large to produce significant damage at the 
center of the plate. At the end of five integration steps, a stress 
analysis was performed. The elements which suffered in-plane and 
interlaminar damage at this time are shown in figures 14 and 15 
respectively. The grid on which the damage is shown encompasses 
all elements which were selected for stress analysis. 

The ply damage shown in figure 14 consists primarily of trans- 
verse (matrix) cracking in the bottom 45® ply. Six of the elements 
also experience matrix cracking within the next interior ply (0®). 
Thus, the region shown in figure 14 has sustained considerable damage 
in terms of surface area but little through the thickness. 

The interlaminar damage shown in figure 15 is more interesting. 
Each of the four elements depicted has experienced delamination at 
the five innermost interfaces. This damage is quite extensive 
and indicates that the transverse shear resultants are extremely 
high in this region and at this time. In practical terms, the 
material in this region must be considered fully degraded. 

The most interesting feature of the interlaminar damage is 
its presence. Had one simply applied a static point load at the 
center of the plate, no delamination would have occurred. However, 
at the very early time at which the stresses were computed, extreme- 
ly sharp bending gradients are present. This is due to the highly 
localized displacement response at this early time, as was shown at 
later times in figures 8 and 9. 

In figure 16, the elements with accumulated damage after 2.5 
- 5 

X 10 seconds are depicted. Comparing figures 14 and 16, the re- 
gion of damaged elements can be seen to be growing. The majority 
of the damaged elements in figure 16 have one or two back face plies 
damaged. Two of the elements shown have suffered top face damage. 
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These elements are shown in figure 17. Obviously at 2.5 x 10 
ceconds, the ply damage has become extensive, with elements fail- 
ing at both outer surfaces. 

-5 

It is interesting that at 2.5 x 10 seconds no additional 
interlaminar damage is present. The reason for this stems from 
the method used for computing interlaminar shear stresses. These 
stresses are computed at ply interfaces only. Since the remaining 
interface was between 0® and 45® plies, the neutral surface will 
not be at the ply interface; the maximum shear stress will also 
not be at the ply interface. Thus, the analysis does not predict 
the maximum shear stress in this case and the lower stress at the 
interface may not cause damage. 

_K 

At 3.75 X 10 seconds, the transverse shear resultants are 
sufficient to produce large interface stresses and additional inter- 
laminar damage is predicted. The new delamination is within the 
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elements which had sustained interlaminar failures at 1.25 x 10 
seconds. Thus, two of the elements shown in figure 15 have only 
one remaining intact p y interface. The remaining interface in 
these elements is between the bottom 0® and 45® plies. 

In figures 18, 19 and 20, the accumulation of ply damage dur- 
ing the remainder of the analysis is shown. No additional inter- 
laminar damage occurred. 

These figures show a steady increase in the region of impact 
induced ply damage. The extent of the damage through the thickness 
at the end of the analysis is depicted in figures 20, 21 and 22. 
Figure 20 indicates the total planar area of the plate damage. Fig- 
ure 21 shows which elements in this region have more than one dam- 
aged ply and figure 22 shows elements which have sustained fiber 
damage. 

The fiber breakage shown in figure 22 is probably the most 

critical. In each of the elements, the first occurrence of fiber 
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damage was at 5 x 10 seconds. Both elements sustained failures 
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in the bottom -45® ply at this time. At 6.25 x 10 seconds, ad- 
ditional damage in the bottom 90® and 0® plies occurred. At 7.5 
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X 10 seconds, each ply in these elements had sustained damage and, 
hence, the solution terminated. There was no material left to carry 
any additional loading. 

The damage accumulated during this solution was very large in 
its extent, both in the planar dimensions of the plate and through 
the thickness. The problem with these results was that the experi- 
ments carried out at NASA Langley did not indicate damage as exten- 
sive as the computer analysis did. The experimental work generated 
interlaminar separation and transverse ply failure as the two pri- 
mary damage modes, but apparently did not produce the extensive fiber 
damage. One reason for this could have been related to the placement 
of the impacting mass. Since the entire mass was lumped at one node 
in the analysis, the response of the plate may have been overestimated. 
It was determined, therefore, that an additional analysis would be 
made with the impact mass distributed over five nodes rather than 
lumped at one node. 

Stress Solution, Mass Distributed at Five Nodes 

The impact mass distribution used in this third analysis was 
previously depicted in figure 4. This distribution was chosen in 
an effort to more closely simulate the impact of a sphere on a plate. 
The diameter encompassed in the distributed mass arrangement corres- 
ponds to 20% of the diameter of the sphere used in the experiments 
at NASA Langley. 

The displacement response predicted in this analysis was very 
similar to the previous analysis. The shapes of the curves were 
nearly identical to those in figures 8 and 9. One surprising dif- 
ference was present, however. The solution with the distributed 
load produced larger peak displacement values. At 2.5 x 10 ‘ sec- 
onds, the difference amounted to 3%. When the time increased to 
-5 

5 X 10 seconds, the difference decreased to 1%. 
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This result was surprising since it was anticipated that the 
spreading o£ the impact mass would reduce the displacements. The 
reason for the increased displacements is the initial velocity solu*- 
tion procedure used in the analysis. Since all nodes where the im- 
pact mass is applied are given an initial velocity, more of the plate 
begins the analysis with non-zero velocities. When only one node is 
impacted, the surrounding nodes lag behind in terms of velocity. 

When the distributed impact mass is used, these surrounding nodes 
are excited at the same velocity. Thus, in the solution where the 
mass is lumped at one node, the velocity lag of the surrounding nodes 
restrains the displacements. This effect is apparently temporary in 
the solution procedure as the difference decreases rapidly with time. 

The distributions of damage in this solution were somewhat dif- 
ferent from those in the previous solution. The extent of the dam- 
age was nearly the same, however, and this solution terminated in a 
fashion similar to the previous case. The increased displacements 
caused this solution to terminate earlier than the lumped mass solu- 
tion. The procedures operated for 20 time steps before the damage 
was too extensive, while the previous solution proceeded for 30 steps. 

At the first stress calculation, the predicted ply and interlami- 
nar damage is depicted in figures 23 and 24 respectively. Comparing 
figures 23 and 14, it can be seen that the extent of the damage is 
very similar for the distributed mass solution and the lumped mass 
solution. The type of damage is also similar in that the majority 
of damaged elements have suffered back ply matrix tensile failures. 

Comparing the interlaminar damage predicted in the two solutions 
(figures 15 and 24) demonstrates that the distributed mass solution 
produces considerably more delamination, although the extent of this 
damage through the thickness of the plate was comparable for the 
two solutions. Each of the damaged elements in figure 24 has suf- 
fered delaminations between all plies except the outer 45®/0® in- 
terfaces. As in the previous solution, this extensive delamination 
must be considered as total failure in these elements. The solu- 
tion procedure continued, however, since most of the individual plies 
remained intact. 
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The damage growth during the remainder of the analysis is de- 
picted in figures 25 through 27. With each stress calculation, the 
number of damaged elements grows, as does the extent of the dam- 
age within the elements. Comparing figures 27 and 19 indicates 
that the number of damaged elements at 5 x 10 seconds is identical 
for the two solutions. The types of damage within the elements are 
very different, however. The damage in the central four elements 
(fig. 28) includes considerable fiber breakage. The solution ter- 
minated because two of these central elements had sustained fiber 
damage in all plies except the upper 0® ply. This ply had suffered 
matrix damage, hov^/ever. Thus, no material remained for carrying 
the load. 

The delamination predicted in the distributed mass solution 

never progressed beyond the twelve elements shown in figure 24. 

The amount of delamination within these elements did increase, how- 
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ever. At 2.5 x 10 seconds, the four central elements (fig. 29) 

had no remaining ply interfaces. During the subsequent stress cal- 

-5 

culation at 3.75 x 10 seconds, four additional elements experi- 
enced increased interlaminar damage. These elements are shown in 
figure 30. The last remaining interface in the elements was the 
bottom 0®/45® ply junction. At the end of the analysis, no addi- 
tional delamination had occurred. 

After reviewing the results of the second stress analysis, it 
was decided that a third should be performed. The decision was made 
in an effort to determine the damage subsequent to the point at 
which the two stress solutions had terminated. 

Stress Solution with Modified Finite Element Model 


In order to continue the solution process beyond the point at 
which an element has no remaining plies, two approaches were con- 
sidered. The first approach involved modifying the computer code 
to simply ignore the element after the damage occurred. The time 
involved in adopting this method was not available, however, and 
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a simple approach was needed. The approach taken involved simply 
removing the four central elements from the finite element model 
before beginning the analysis. It was felt that this approach would 
eliminate the most critical elements and allow the solution to pro- 
ceed to conclusion. Since these four elements had already been 
shown to sustain damage throughout the laminate, removing them would 
simulate the response after they had failed. 

The finite element model used for this analysis was previously 
shown in figure 3 and the impact mass distribution demonstrated in 
figure 4 . 

The solution with this modified finite element model progressed 
through 30 time steps, as had the original stress solution. After 
30 steps, the solution again terminated due to the lack of material 
left in an element. The progression of damage growth for this so- 
lution is shown in figures 31 through 37. The only major differ- 
ence between this solution and the previous ones relates to the 
location of the element which caused the solution to terminate. 

In the previous analyses, the element which failed totally had 
always been one of the four central elements. In the modified model 
solution, these elements were removed. The element which caused the 
execution to terminate in the modified model solution was simply an 
adjacent element. Hence, the effort to continue the analysis be- 
yond 7,5 X lO”^ seconds was in vain. 
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IV. DISCUSSION AND RECOMMENDATIONS 


The computerized analysis methodology has been utilized to 
predict the response r including damage initiation and propagation, 
for a specific low velocity, low mass impact event. The analysis 
has been shown to be an effective tool in predicting the response 
of the plate structure. The damaged nodes and locations predicted 
were consistent and correspond reasonably well with experimental 
work as described to MSC by NASA Langley. Several significant 
features of the analysis predictions and methodology warrant fur- 
ther discussions. 

In each of the stress solutions performed, significant inter- 
laminar separations occurred at the first stress computation. The 
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time at which the calculations were made was 1.25 x 10 seconds. 

The elements which suffered interlaminar separation at this time 
were the only elements affected by delamination. These two features 
identify two significant elements of the stress and dw.'nage predic- 
tions . 

The fact that the delamination occurred at the first stress solu- 
tion demonstrates that interlaminar failures initiate very early in 
the impact event. The delamination is a function of the transverse 
shear forces, which are highest very early in the impact event be- 
cause of the localized nature of the displacement response at that 
time. Since the significantly non-zero displacements are restricted 
to a small region at the center of the plate, severe moment gradients 
also exist at the plate center. The high moment gradients produce 
the large transverse shear forces. As the impact event progresses, 
the displacement fields spread, and the moment gradients decrease 
and the propensity for interlaminar damage initiation decreases. 

As delamination occurs, shear deformations increase due to 
singular shear stresses at the crack tip. The singular shear 
stresses tend to promote propagation of the crack. These effects 
cannot be modeled, however, since shear deformation is not in- 
cluded in the current analysis. 
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This is an area where further work could significantly improve the 
usefulness of the CLIP code. The propagation ot interlaminar dam- 
age might be handled in one of two ways. 

One approach would involve the inclusion of shear deformation 
in the displacement solution. With the inclusion of shear defor- 
mation, damaged elements could be given a reduced stiffness to 
promote increased shear stresses in adjacent elements equal to 
the average increase produced by the shear stress singularity. 

This would then promote delamination propagation. 

The other approach would involve the use of an analysis to be 
used subsequent to the impact analysis. This method would use the 
delamination predicted during the impact analysis. The material 
surrounding the delamination would then be subjected to the shear and 
moment distributions predicted in the impact analysis after the 
instant at which the damage occurred. The areas where the delamina- 
tion connected with intact material could then be analyzed utilizing 
fracture mechanics to predict the delamination growth. Noting that 
the effects of shear deformation are small in undamaged materials, 
this damage growth need not be included in the impact analysis if 
the affected area is small. 

Another area which warrants discussion involves the extent of 
the ply damage predicted in the impact analyses made. In all of the 
stress solutions made, the damage eventually propagated completely 
through the thickness of the central finite elements. This led to 
an immediate problem in that the solution terminated at this point. 
Modifications to the CLIP code to allow the solution to proceed 
beyond this point would be most helpful. The solution could be modi- 
fied to continue until the damaged region reached an area equivalent 
to the impacting mass. At this point one would assume that the mass 
had penetrated the plate. 

Because the solutions terminated, the full extent of the damage 
which would have been induced in the plate was never determined. It 
was apparent, however, that the predictions of damage exceeded the 
damage measurec in the experimental work performed at NASA Langley. 
This was primarily related to the damage through the thickness 
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predicted at the plate center and this would tend to indicate that 
the energy impacted to the system in the analysis was too great. 

The excess energy in the analysis is probably a function of two 
effects not modeled. The first of these, and probably the most 
important, is related to energy which should be lost due to local 
surface crushing at the point of impact. It is known from Hertzian 
contact analyses that the local contact stresses can be quite large 
and could cause significant permanent local deformations. These 
deformations would absorb a significant portion of the impact energy. 
Modifications to account for this effect would improve the analytical 
predictions. These modifications could be effected utilizing the 
contact forces prediction already in the CLIP code coupled with a 
Hertzian contact analysis and a local, non-linear material model. 

The other contributor to the energy loss is material damping. 

The effects of damping very early in the impact event are probably 
small but at later times they must surely become significant. The 
CLIP code has provisions for damping but this feature was not used 
due to a lack of data. 

Another area of discussion relates to the effects of utilizing 
a distributed impact mass in the analysis. As was demonstrated, 
distributing the impact mass produced a slightly increased dis- 
placement field when the opposite should have been true. The effect 
was seen to be short-lived and may be insignificant. It is discon- 
certing, however, and modifications could be made to correct the 
situation. These modifications would involve applying the mass and 
velocity initially at one node and spreading the mass as a function 
of the contact force as the solution progressed. This would more 
closely simulate the actual impact event. 

A final area of consideration relates to the effects of a 
static prestress on the impact analysis. The procedure utilized in 
the current analysis is simply a superposition of the static and 
dynamic displacement fields. This is applicable if in-plane static 
loads are tensile or if shear or compressive loadings are small 
enough that buckling is not a consideration. If buckling is a real 
possibility then the inclusion of a stability analysis would be 
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required . 

A etabllity analysie in the finite element procedure would 
require considerable modifications to the code and substantially 
decrease the efficiency of the analysis. A stability analysis 
requires the formuxation of a modified stiffness matrix and an 
eigenvalue extraction. The modified stiffness matrix .Is based on 
the strain field present when the buckling analysis is performed. 
TI .2 eigenvalue problem is not compatible with the direct integra- 
tion of the equations of motion as currently used for the impact 
displacement response. Thus, a new analysis procedure would be 
needed. Additionally, the eigenvalue problem would have to be 
solved at each time step since the buckling loads are a function of 
the strain field present. This would effectively double the solu- 
tion time. Thus, while the inclusion of a stability analysis would 
be possible, it would not be practical. 

In summation, while the analysis is functioning well and pro- 
vides significant insight into the phenomenon of low velocity im- 
pact, further work would be most helpful. 

The areas of these future efforts should include: 

1) Inclusion of a disbond growth model either 
thx'ough shear deformation or fracture mechanics; 

2) Modification to allow continued program execution 
after local complete element failure; 

3) Inclusion of energy loss mechanisms due to non- 
linear contact effects; and 

4) More realistic modeling of the impact mass 
distribution. 

The addition of these capabilities to the analysis procedure 
would provide better modeling of the impact event. This would al- 
low for increased confidence in the analysis and facilitate its use 
in evaluating the relative merits of various laminates when sub- 
jected to low velocity impact. 
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A8 optimum laminates are determined from impact related cri- 
teria, further evaluations could then be performed with regard to 
residual static strength, subsequent impact events, and fatigue 
The static and fatigue evaluations could be made using a 
modification of a code developed for fatigue of notched composite 
laminates (reference 18) while subsequent impact events might be 
evaluated using a modification of the code developed here. 
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Table 1. Failure Criteria 
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Table 1 (cont'd.). Failure Criteria 


Tensile Interlaminar Mode (^33 ^ 0.0) 
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Compressive Interlaminar Mode < 0.0) 
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where : 

axial tensile strength 
= axial compressive strength 
= transverse tensile strength 
0 ^ = transverse compressive strength 
= axial shear strength 
- transverse shear strength 


Table 2. T300/5208 Properties 


153 GPa 
10.9 GPa 

5.6 GPa 
0.30 

1.55 g/cc 
0.70 

689.5 MPa 

758.5 MPa 

27.6 MPa 

96.5 MPa 
62.1 MPa 
62.1 MPa 
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Table 3 


Impact Parameters 


IMPACT MASS = 16.45 g 
IMPACT VELOCITY =9.4 m/s 
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Table 4. Comparison Groupings of Predicted Damage 


Figures 14 - 

22 

Predicted damage in solution with mass lumped 
at one node 

Figures 23 - 

30 

Predicted damage in solution with mass dis~ 
tributed over five nodes 

Figures 31 - 

37 

Predicted damage in solution with modified 
finite element model (four central elements 
removed ) 

ONE NODE IMPACT 

Figures 


Comparison 

14, 16 - 20 


Figures show the extent of elements with 
predicted ply damage as a function of time 

14, 15 


A comparison of the prediction of elements 
with ply damage vs. elements with interlami- 
nar damage at T = 1.25 x 10'"^ sec 

16, 17 


A comparison of the number of elements with 
predicted top ply damage vs. the number of 
elements with ply damage anywhere through the 
laminate thickness at T = 2.5 x 10“^ sec 

20 - 22 


A description of the extent of ply damage 
through the thickness of the laminate 

FIVE NODE IMPACT 

Figures 


Comparison 

23, 25 - 27 


Figures show the extent of elements with 
predicted ply damage as a function of time 

23, 24 


A comparison of the prediction of elements 
with ply damage vs. elements with interlami- 
nar damage at T = 1.25 x 10"^ sec 

27 - 30 


Figures depict the extent of damage in the 
laminate including ply damage, fiber damage 
and interlaminar damage 
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Table 4 (continued) . Comparison Groupings of 
Predicted Damage 


MODIFIED FINITE ELEMENT MODEL 


Figures 


Comparison 


31, 33, 34 - 37 Figures show the extent of elements with pre 

dieted ply damage as a function of time 


36, 37 A comparison of elements with predicted ply 

damage vs. elements with interlaminar damage 
at T = 1.25 X 10“^ sec 


Figures 
15, 24, 32 


19, 27, 35 


20, 27, 37 


20, 28 


COMPARISONS BETWEEN ANALYSES 

Comparison 

A comparison of the number of elements which 
have sustained interlaminar damage in each 
of the three solutions 

A comparison of the number of elements which 
have sustained ply damage at T = 5.0 x 10“^ 
sec in each of the three solutions 

A comparison of the number of elements with 
predicted ply damage at the end of each solu- 
tion 

A comparison of the number of elements which 
have sustained fiber damage in the solution 
with one impacted node vs. the solution with 
five impacted nodes 
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15.24 


X 



15.24 



Figure 3. 
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i 

c) Modified 



Figure 4. Impactor Mass Distributions 


















Y 

max 


Figure 10. Moment Resultants Through the Plate Center Along 
the Shorter Axis, Stress Solution, Mass Lumped at 
One Node . 
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Figure 11 . Moment Resultants Through the Plate Center Along 
the Longer Axis, Stress Solution, Mass Lumped at 
One Node . 
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Figure 14 . 


Elements With Ply Damage, T = 1.25 x 10 
sec., Stress Solution, Mass Lumped at One Node 
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Figure 19. Elements with Ply Damage, T = 5.0 x 10 sec.. 
Stress Solution, Mass Lumped at One Node 
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Figure 20. 
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Elements with Piy Damage, T = 7.5 x 10 sec. 
Stress Solution, Mass Lumped at One Node 
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Figure 21. Elements withgMore Than One Failed Ply 
T = 7.5 X 10 sec., Stress Solution, 
Lumped at One Node 
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Figure 23. Elements with Ply Damage, T = 1.25 x 10 

sec., Stress Solution, Mass Distributed at 5 
Modes 
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Figure 24. Elements with Interlaminar Damage, T = 1.25 x 10 

sec.. Stress Solution, Mass Distributed at 5 Nodes 
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Figure 26. Elements with Ply Damage, T = 3.75 x 10 sec., 
Stress Solution, Mass Distributed at 5 Nodes 
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Figure 28. Elements with Fiber Damage, T = 5.0 x 10“ s 
Stress Solution, Mass Distributed at 5 Nodes 


ORIGINAL PAGE IS 

OF POOR QUALITY 




Figure 29. Elements with No Remaining Intact Ply Interfaces, 

T = 2.5 X 10”5 sec., Stress Solution, Mass Distributed 
at 5 Nodes 
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Figure Elements with Ply Damage, T = 3.75 x 10 

sec., Stress Solution, Modified Model 
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Figure 35. Elements with Ply Damage, T » 5.0 x 10 
sec., Stress Solution, Modified Model 
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APPENDIX A 


STRESS, MOMENT, AND TRANSVERSE SHEAR RESULTANTS 

In order to perform the composite stress analysis required for 
failure prediction, it is necessary to determine stress, moment, 
and transverse shear resultants. The thin shell finite element 
taken from SAP IV has the capability for stress and moment result- 
ants but lacks the ability to compute transverse shear result? its. 

It was necessary, therefore, to add the capability for this trans- 
verse shear prediction. 

Since the primary focus of the current contract involved ma- 
terial response and not the development of general purpose finite 
element routines, the .implest approach available for computing 
transverse shear resultants was adopted. In the course of this 
effort, it was determined that since the restrictions inherent in 
the simple transverse shear computations effectively remove the ne- 
cessity for the complex stress and moment resultant computation, 
the SAP IV routines which perform these analyses were not required. 
Therefore, the entire stress recovery procedure is based on the 
simplifying element orientations required for transverse shear 
as described below. 

The primary restriction invoked by the simple transverse shear 
resultant computation involves the elemental geometry. No attempt 
was made to determine transformations of forces from the global 
coordinate system to an arbitrary elemental coordinate system. 
Therefore, it is required that the local elemental coordinates coin- 
cide with the global coordinates. Additionally, no provisions were 
made for elemental geometric irregularities. This requires that all 
elements be rectangular. The second restriction applies only to the 
stress recovery procedures. The displacement response is computed 
correctly for elements which are not rectangular if the element co- 
ordinate system is defined coincident with the global coordinate 
system. This restriction is required so that the laminate prop- 
erties are utilized correctly. 

In figures A-1 and A-2, the elemental nodal forces and moments 
are depicted. The forces and moments are shown in a global positive 
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scene. The element is required to be rectangular and oriented with 
the global coordinate system as mentioned previously and^ thus, the 
relations for stress and moment resultants can be written directly 
as given in tables A-1 and A-2. It can easily be seen that the 
relations in tables A-1 and A-2 are invalid if either of the two 
restrictions mentioned are violated. 

The computation of transverse shear stress resultants is slight- 
ly more complicated. In figure A-3 the moments and shear resultants 
required for equilibrium are depicted. In the analysis used, the 
transverse shear forces are computed from moment equilibrium since 
shear deformation is not included in the displacement solution. 

The transverse shear resultants can be determined quite easily 
from equilibrium considerations as: 


3M. 


3M 


Q. 


3x 




ay 


Q.. = - 


3M 3M 

— iai z 

3x 3y 


(A.l) 


If it is assumed that all quantities vary linearly within the 
finite element, equations [A.lj can be rewritten as: 





AM. 


2 ^ 


Ay 




Ax 


AM 

^ . 

Ay 


(A. 2) 


Each of the quotients in equations [A. 2] can then be written in 
term:, of the nodal moments of figure A-2. When this is done, a dif- 
ficulty is immediately encountered. 
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The various terms of equation [A. 2] are found to be: 


Ax 

AM 




AM 

xy 

“Sx 


AM 


xz 


Ay 


- [mJ + mJ + + mJ] 

imJ; + + M^] 

[m’- + m^ + m^ + m^] 

X X X X 

- (M^ 4 M^ + M^ + M^] 

y y y y 


1 

AxAy 



AxAy 


_1 

AxAy 


1 

AxAy 


(A. 3) 


It is immediately apparent that in terms of the nodal moments, 
the two terms required for Q are identical, as are the two terms 
in the expression. This result is demonstrating that, in terms 
of the nodal moments, it is not possible to separate the contribu- 
tions of the bending and twisting moments. This is easily confirmed 
by considering a one-dimensional case. The shear force computed is 
twice the proper value if the full expressions of equations [A. 2] 
are utilized. 

The resolution of this problem involves simply dropping the 
twisting components from equation [A. 2] and using the expressions; 


AM 


Q.. = 


x 


Ax 




AM 

- __z 

Ay 


(A. 4) 
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It must be remembered that the contributions of the twisting 
moments have not been neglected in equations [A. 4]. These contribu- 
tions are included implicitly in equations [A. 4] through the compu- 
tations of the two quotients involved. 


Table A-1. Stress Resultants 


N = [ (F^ + F^) 

X X X 

N = [(F^ + F^) 

y y y 

N = [(F^ + F^) 

xy X X 

= [(F^ + fJ^) 


(F^ + F^)] 
X X ^ 


(F^ + F^) ] 

y y 

(F^ + F^) ) 
X X 

(F^ + F^) ] 

y y 


1 _i 

2 Ay 


1 _1 
2 Ax 


1 J, 

2 Ax 


1 _1 
2 Ay 
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Table A-2 


Moment Resultants 


M 


X 




(<mJ + mJ) 
[(m‘ + mJ^) 
[(m!) + 

X X 

[(M^ + M^) 

y y 


(M^ + mJ^) 1 
j y 

(M^ + M^) 1 
X X 

<“x * <>> 

(M^ M^) ] 

j y 


1 1 

1 Ay 


1 _l 

2 Ax 


1 _1 
2 Ay 


1 _1 

2 Ax 
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•H X 



A”2. El6inGnt3l Nodsl Morncnts 
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APPENDIX B 


INTERLAMINAR SHEAR STRESSES 


The analysis used in the CLIP code for computing interlaminar 
shear stresses is based upon bending equilibrium. In Appendix A, 
the derivation of the transverse shear stress resultants is given. 
These resultants are then used to compute transverse shear stresses 
at the ply interfaces within the laminate. The method used for com- 
puting the interlaminar shear stresses is based on an extension of 
classical methods used for bending shear stresses in homogeneous 
beams . 

In the analysis, it is necessary to compute interlaminar shear 
for two different cases. Interlaminar stresses must be computed in 
elements which have sustained damage as well as those which have 
not. The methods used in both cases are described here. 

INTERLAMINAR SHEAR IN INTACT ELEMENTS 

In order to compute interlaminar shear stresses in undamaged 
elements, the transverse shear stress resultants are converted to 
moments by multiplying by the appropriate shell element dimensions 
(see eqn. A. 4). These two moments then represent the change in 
moment across the element in both the X and Y coordinates. This is 
demonstrated for a one-dimensional case in figure B-1. It is readily 
apparent in figure B-1 that the shear stresses represent the balanc- 
ing force required for moment equilibrium. 

The moment differentials derived from the transverse shear 
resultants are then applied to the laminate model. This produces 
three in-plane stresses within each ply. These ply stresses are 
then converted to forces and summed through the thickness in both 
the X and Y coordinate directions. In figure B-2, the three in- 
plane stresses and resulting interlaminar shear stresses are shown 
for an outer ply of the laminate and an applied moment differential. 
Equilibrium in the X-direction yields forces 
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0^ Ay Lz + a Ax Az + o Ax Ay ■ 0 
X xy X z 


(B.l) 


and in the Y-direction 


o„ Ax Az Ay Az + a Ax Ay » 0. 

y y* 

Thus, the Interlaminar shear stresses are simply 

-(o Ay Az + 0 Ax Az) 

X ^y 

CJ 8 ' '■ I ■!. ■ .1 III 

Ax Ay 


a 


yz 


(a Ax Az + a Ay Az) 
y 

Ax Ay 


(B.2) 


(B.3) 


To compute interlaminar shear stresses on interior ply inter- 
faces, it is only required that the X and Y forces represented by 
the and terms respectively in equation B.l be summed through 
the laminate to the interface in question. 


INTERLAMINAR SHEAR IN DAMAGED ELEMENTS 

In order to compute interlaminar shear stresses in an ele- 
ment which has sustained damage prior to the stress calculation, 
it is first necessary to determine how the laminate responds to 
loading when damage is present. In figure B-3 the bending 
stresses through a localized delamination are shown. Since the 
delamination is smal 1 , the material away from the delamination 
effectively forces the curvature and strain field to remain un- 
changed. The material above and below the delamination bend as 
two independent laminates while the constraint of the adjacent, 
undamaged material adds opposing membrane forces. Thus, the net 
effect produces no change in the bending stress field. This is 
the basis for not reformulating the elemental stiffness for 
interlaminar delaminations. 
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This model cannot account for a moment gradient, however. In 
order to evaluate the effects of a mo'fent gradient, it is required 
that the two sublaminates behave entirely independently. The op- 
posing membrane forces must not be present in the shear analysis. 

If a moment gradient is applied to this type of model, then the mem- 
brane forces must also exhibit a gradient. This is not possible 
since it would require a shear force transfer across the delamination. 
This is a drawback related to the lack of shear deformation in the 
analysis . 

To model the moment gradient, it is necessary to compute the 
bending stiffnesses of the sublaminates independently. The bending 
stiffnesses are then added to give the stiffnesses of the assemblage 
of partial laminates. This, however, leads to bending/extensional 
coupling at both the laminate and sublaminate levels. It is required 
therefore, that the [B] matrices be accounted for. This can 
easily be accomplished by noting that the membrane forces must be 
zero for the applied transverse shear forces (moment gradients) . 

The force/moment - strain/curvature relations for a general 
laminate are 


(N) = + [bHk) 

(M) " lB]{e°) + [dHk) 


(B.4) 


where the matrices have the usual connotations. Remembering that 


(n) = 0 


(B.5) 


the mid-plane strains due to a bending load are 

= -[A]"^ [B] {<} . (B.6) 

These strains actually represent a neutral axis shift from the 
mid-surface of the unsymmetric laminate to the proper neutral surface 
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Using the strains of equation (B.6) » an effective bending 
stiffness matrix can be determined. 


lBj{e°} + lD]{i<} 

(B.7) 

-IB] IA]“^ IB] {<} + ID] {«} 

(B.8) 

(ID] - IB] IA]‘^ IB]) ( k ) 

(B.9) 

ID*] (k) 

(B.IO) 


The bending stiffness matrix in (B.IO) relates moments and curva- 
tures about the neutral surfaces rather than the mid-surface. 

The bending stiffness matrix for the total laminate is then 
given as 

^ I 

ID*] = E ID*] , N - number of sublaminates. (B.ll) 

1=1 

With the bending stiffness of the assemblage of sublaminates given 
in eqn. (3.11), a curvature vector can easily be found for the ap- 
plied moment differentials. Using the curvature for the whole lami- 
nate and eqn. (B.6), the ply stresses for each sublaminate can be 
determined and summation of forces depicted in figure B-2 and equa- 
tions B.l through B.3 carried out. Thus, the interlaminar stresses 
in damaged elements are predicted. 

When an element contains both ply damage and interlaminar damage, 
the same method is used. In figure B-4 the sublaminates for a lami- 
nate which has sustained general damage are shown. For ply damage, 
the entire ply is removed from interlaminar shear calculations. 
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Figure B-l. 


One-Dimensional Shear/Moment Relationship 
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Figure B-4. Typical Failed Laminate 





APPENDIX C 


COMPOSITE LAMINATE IMPACT PROGRAM (CLIP) 


The computer program developed for predicting damage initiation 
and growth during low velocity impact consists of a transient dynarn- 
ics finite element code coupled with composite stress and failure 
analysis procedures. In conjunction with the failure analysis rou> 
tines ifi the capability to incorporate predicted damage into the 
transient analysis. Thus, the damage predicted during any time step 
is incorporated into the dynamic solution at future time steps. 

This coupling of damage and dynamic response is the heart of the com- 
puterized procedure. 

The transient dynamic finite element procedures used were taken 
from the SAP IV code (ref. 14) . The portions of the SAP IV code 
used include the thin shell element routines and the time integra- 
tion routines. The composite stress and failure analysis routines 
were developed specifically for use in the CLIP code. 

In figure C-1, a simplified flow chart of the analysis performed 
by the CLIP code is shown. The branches and loops shown represent 
the procedures required for incorporating the effects of predicted 
damage into the analysis and removing the impact mass after the 
contact force becomes tensile. The various computer program seg- 
ments which perform the analysis are listed in table C-1 along with 
a brief description of their functions. 

PROGRAM OPERATION 

The CLIP code has been written with several different solution 
options. The different options offer the user considerable flexi- 
bility in using the code. The options available include data check 
mode, static analysis alone, static analysis to generate a pre- 
stress condition for a dynamic analysis, dynamic analysis without 
stress calculation and dynamic analysis with stress calculation. 

The data check mode of operation is useful in insuring that 
all input data parameters are correct. This option is essential 
in any large program. 
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The static solution mode has, basically, tv;o uses. First, a 
static test case is an easy method to prove the validity of a finite 
element model, it is extremely difficult to verify the results of 
a dynamic analysis directly. The second use of the static analysis 
is to generate a pre-stress condition at the beginning of the dynamic 
analysis. This allows the static stresses to be included directly 
with the dynamically induced stresses. 

The inclusion of the static solution is performed by direct 
superposition. The static displacement vector is saved at the end 
of each time step. In this way the static displacements do not act 
as an initial condition in the dynamic response. It must be remem- 
bered that this is a direct superposition of results. There is no 
coupling of the static and dynamic responses. 

The dynamic solution can be performed with or without stress 
calculations. The dynamic analysis made without stress calculations 
is considerably less costly and time consuming than a similar analy- 
sis with a complete stress analysis. The solution cost reductions 
are a result of two factors. First, the stress analysis procedures 
are by-passed in the analysis and secondly, there are no damage cal- 
culations which would require reformulation of elemental and global 
stiffness matrices and subsequent decomposition. Because of this, 
the dynamic solution without stress analysis can be used to determine 
the length of time required to characterize the impact event at a 
moderate cost. This mode of analysis can also be used to verify 
the integration time step selected. 

When using the dynamic solution without stress analysis, it 
should be remembered that the response of the plate will be different 
when stress calculations are included. The amount of this difference 
will be a function of the amount of damage sustained in the plate 
when stress calculations are included. 

When damage predicted in an element is ply damage {either matrix 
or fiber) , the elemental stiffness matrix is modified. These modi- 
fications are then included in the global stiffness matrix for in- 
corporation in subsequent time steps. Thus, the plate stiffness is 
changed forcing a plate response change. 
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The dynamic solution mode which includes stress analysis is the 
primary mode of operation of the CtIP code. This mode is used for 
simulating the actual response of a composite plate subjected to low 
velocity impact. 

The use of the analysis in this mode requires some thought on 
the part of the user. The code is written such that the user selects 
how often the stress calculations are performed. If the stress cal- 
culations are not performed often enough , it is possible to miss 
considerable damage. This is especially true at the beginning of 
the analysis when very localized deformations are the predominant 
response of the plate. 

In an effort to achieve computational economy, the CUP code 
performs stress analysis only for elements which the user selects. 

This option saves considerable computet time since many of the 
finite elements in the model will typically never sustain damage. 

If, however, an element, which should fail under the dynamic load- 
ing, is not selected for stress analysis, this damage is lost. It 
is never computed and, hence, is not included in the dynamic response. 
Thus, it is important to select the proper group of elements for 
stress calculations. Insight into the elements which should be in- 
cluded can be obtained by examining the plate response without stress 
calculations. 

PROGRAM OUTPUT 


The output of the CLIP code can be divided into two distinct 
entities. First, the program echos the input data such that a 
record of this data is generated, and secondly, the results of the 
analysis are printed. In each of these two areas, the user has 
considerable flexibility in determining the extent of the output 
generated. 

In terms of the input data echo, the user can suppress printing 
of the nodal points or the elemental data or both. Other input data 
including the impact parameters, lamina materials and laminate con- 
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figuratior as well as data generated by the analysis code with re- 
spect to problem size and laminate elastic constants are always 
printed. 

In terms of the solution output, the user selei'ts the frequency 
of output for the various data generated as well as the nodal points 
(displacements, velocities, accelerations) and elements (element 
forces, laminate coordinate stresses, layer coordinate stress ss, 
failure calculations, failure locations) where data are to be printed. 

For the nodal output data, the user can select different print- 
out frequencies for the displacements, velocities and accelerations 
as long as the velocity and acceleration print frequencies are even 
multiples of the displacement frequency. Thus, the displacements 
may be printed every five time steps while the velocities print 
every fiftj steps and the accelerations every 20 steps. The print 
frequencie;* are left to the discretion c£ the user. 

The elemental data are printed in a similar fashion where each 
printout frequency must be an even multiple of the stress calculation 
frequency. As an example, the user must not select a stress calcula- 
tion frequency of every ten steps while requesting a layer coordinate 
stress printing every three time steps. The result of this particu- 
lar arrangement would produce layer stress printing only every thirty 
time steps. There is no requirement that the individual data print 
frequencies be multiples of each other, however. 

The elemental data which can be printed consist of elemental 
forces, ply stresses in two coordinate systems, calculations of the 
failure analyses, and locations of damage. The elemental forces 
consist of stress, moment and transverse shear resultants and top 
surface stress corresponding to the impactor/plate contact force 
distributed over all elements adjoining the impacted nodes. 

Layer stress can be printed in the laminate coordinate system 
(global X-Y) or the local ply coordinate systems, or both. Stresses 
are printed for each ply in the laminate under study as well as each 


98 


ply interface including the top and bottom free surfaces. Thus, 
there is always one more interface than there are plies. 

The calculations of the failure criteria, if selected for 
printing, produce data for each ply and interface. The data 
printed are simply the sum of the terms in the various failure 
criteria modes (see table 1). These data can be used to iden- 
tify locations and modes of damage within the elemental laminate 
models. 

If the details of these calculations are not required, the 
user can select that failure location printing he activated. 

This output details the ply or interface which has failed within 
a damaged element. For the case of ply damage, information is 
also printed indicating whether matrix or fiber failure has 
occurred. It is possible to specify both failure calculation 
printing and failure location printing. To do so, however, is 
som ,‘^wha t redundant . 
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CLIP USERS GUIDE 


Composite Laminate Impact Program 
I. PROGRAM CONTROL DATA 



Columns 


Contents 

1-5 

NT 

Maximum number of time steps 

6-15 

DT 

Time step 

16-25 

ALFA 

Mass proportional damping 
coefficient 

26-35 

BETA 

Stiffness proportional damping 
coefficient 


100 



Card 4 


Print Control Data 


1015 


Columns 



Contents 

1-5 

KEY(l) 


Node print code, 
print node data 

6-19 

KEY(2) 


Element print code^ print 

element data 

11-15 

KEY (3) 


Displacement print code 

16-20 

KEY(4) 


Velocity print code 

21-25 

KEY(5) 


Acceleration print code 

26-30 

KEY(6) 


Element force print code 

30-35 

KEY(7) 


Laminate coordinate stress 
print code 

36-40 

KEY(8) 


Lamina coordinate stress 
print code 

41-45 

KEY (9) 


Failure calculation print 
code 

46-50 

KEYC 0) 


Failure location print code 

NOTES : 

Card 4 




KEY (3) -KEY (10) 

•>0, suppress printing 

»N, print every N time steps 


KEY(6)-KEY(10) 

Codes 

not used if NSTR, Card 2, ^^O 

Card 5 


Print 

Control Data 2015 

Column 



Contents 

1-5 

NDISP 


Number of nodal print groups £100 

6-10 

NSTRP 


Number of element stress 


calculation groups £100 
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N OTES ; 


Card 6 

Colum n 

I- 5 
6-10 

II- 15 
16-20 


Card 7 

Column 

I- 5 

6-10 

II- 15 
16-20 


NOTES : 


Card 5 

If NOISP ■Of print data for all nodes, 
do not include Card 6 

If NSTRP *0, compute stresses for all elements, 
do not include card 7 

Print Control Data 2015 




Contents 

KEYS 

(1,1) 

First node, print group I 

KEYS 

(1,2) 

Last node, print group I 

KEYS 

(J,l) 

First node, print group J 

KEYS 

(J,2) 

Last node, print group J 


Continue through NDISP groups, more than 1 card 
if nec/ii«3ary 


Print Control Data 2015 


KEYS (1,1) 
KEYS (1,2) 
KEYS (J,i) 
KEYS (J,2) 


Contents 

First element stress calculation, 
Group I 

Last element stress calculation, 
Group I 

First element stress calculation, 
Group J 

Last element, print group J 


Continue through NSTRP groups, more than 1 card 
if necessary 


Cards 5,7 


Stress calculations are performed only for element 
identified in these groups. 
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II. IMPACT CONTROL DATA 


Card 1 


Node Data 


15 


Columns Contents 

1-5 NUMIMP Number of impacted nodes 

0 < NUMIMP < 20 


Card 2 


Node Data 

2015 

Columns 



Contents 


1-5 

IMN(I) 


Impacted node 


0 

1 

IMNd-*-!) 


Impacted node 



Continue 

through NUMIMP 

nodes. 


Card 3 


Impact 

Conditions 

2F10.0 

Columns 



Contents 


O 

1 

TOTMAS 


Total impactor mass 


11-20 

TERMV 


Impact velocity 


Card 4 


Impact Mass Distribution 

8F10.0 

Columns 



Contents 


O 

1 

XMFRAC ( I ) 


Impact mass fraction, 

node I 

11-20 

XMFRAC(J) 


Impact mass fraction, 

node J 


Continue 

through NUMIMP 

nodes . 
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IIL; LAMINATE DATA 


Card 1 


Laminate Geometry Control Data 215 

Columns 


Contents 

1-5 

NMAT 

Number of different 



materials < 2 

6-10 

NPLY 

Number of plies £ 25 

Card 2 


Lamina Elastic Constants 5F10.0 


Columns 

I- 10 El (I) 

II- 20 E2(I) 

21-30 G(I) 

31-40 ANU(I) 

41-50 RHO(I) 

Continue through NMAT 
Notes : Card 2 


Contents 

Lamina Axial Modulus En, 
material I 

Lamina Transverse Modulus, 

E 22 # material I 

Lamina Axial Shear Modulus, 

Gi2t material I 

Lamina Major Poissons Ratio 

^^12, material I 

Lamina Mass Density, material I 
rials . 


Major Poissons Ratio, vj[2 * 

Card 3 Lamina Material Strengths 6F10.0 


Columns 

I- 10 STREN(1,I) 

II- 20 STREN(2,I) 


Contents 

Axial tensile strength, material I 
Axial compressive strength, 
material I 
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1 


21-30 

STR£N(3,I) 

Transverse tensile strength, 
material I 

31-40 

STREN(4,I) 

Transverse compressive strength 
material I 

41-50 

STREN(5,I) 

Axial shear strength, material 

51-60 

STREN(6,I) 

Transverse shear strength, 
material I 


Continue 

through NMAT materials. 

Card 4 


Lamina Geometries 2F10.0,I5 

Columns 

Contents 

1-10 

THIK(I) 

Ply thickness 

11-20 

THET(I) 

Ply orientation, degrees 

21-25 

MAT (I) 

Ply material number 


Continue through NPLY plies. 
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IV. NODAL DATA 


Card 1 Nodal Input Data 7I5«3F10.0 

Columns 

I- 5 N 

6-10 ID(N,1) 

II- 15 ID(N,2) 

16-20 ID(N,3) 

21-25 ID(N,4) 

26-30 ID(N,5) 

31-35 ID(N,6) 

36-45 X(N) 

46-55 Y(N) 

56-65 Z(N) 

Continue through NUMNP nodes 
NOTES ; Card 1, 

ID(I,J)=0, Force boundary condition 

= 1, ZeiTO displacement boundary condition 
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Contents 

Nodal point number 
Boundary condition code, 
global X - direction 
Boundary condition code, 
global Y - direction 
Boundary condition code, 
global Z - direction 
Boundary condition code, 
global X - rotation 
Boundary condition code 
global Y - rotation 
Boundary condition code 
global Z- rotation 
X - Coordinate 
Y - Coordinate 
Z - Coordinate 



V. ELEMENT DATA 


Card 1 Eleme>it Input Data 15 

Columns Contents 

1-5 NUMEL Number of elements 


Card 2 Element 


Columns 

I- 5 MM 

6-10 IY(1) 

II- 15 IY(2) 

16-20 IY(3) 

21-25 IY(4) 

26-30 IY(5) 


31-40 PRESSU 


i nput data 6I5,F10.0 

Contents 
Element number 
Element node I 
Element node J 
Element node K 
Element node L 
Element stiffness re-use 
code 

«0, form new stiffness 
=■1, re-use previous stiffness 
Uniforn lateral pressure 
load 


Continue through NUMEL cards. 

NOTES ; Card 2 

IY(1) - IY(4), 

1. Element must be rectangular 

2. Nodes IY(1), IY(2) must lie along global X-Axis 

3. IY(1)-IY(2) Define local positive X-Axis 

4. Element must be defined counter-clockwise 


IY(5), for stiffness matrix re-use. 

Element geometry must be identical to 
previous element 
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PRESSUf Used in static analysis only 

Ignored if ISTAT Section I, Card 2 
VI STATIC NODAL LOADS 


Card 1 


Nodal Load Input 


I5,F10.0 


Columns 
1-5 N 

6-15 TR(1) 

16-25 TR(2) 

26-35 TR(3) 

36-45 TR(4) 

46-55 TR(5) 

NOTES : Card 2 


Contents 
Node number 
X-Axis load 
Y-Axis load 
Z-Axis load 
X-Axis moment 
Y-Axis moment 


Used in static analysis only, but at least one blank 
card must be included. 


N - If N*0 or blank, terminate nodal load input 

Z- Axis moments are not permitted since only flat 

models are allowed and thus no Z-rotational stiffness 
exists 


- END OF INPUT DATA - 
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GENERAL PROGRAM NOTES: 


1. Units for the various input parameters need only be consistent, 
with the exception of ply orientation angles, which must be 
input as degrees. 

2. When inputting nodal boundary condition codes, it is advis- 
able to use zero-displacement wherever possible. These 
serve to reduce program size and running time. Specifically, 
Z-axis rotations should be excluded. X and Y axis displace- 
ments should be excluded unless membrane forces are included 
in static prestress conditions. 

3. Z - coordinate location must be identical for all nodes. 

This restriction, along with che restriction that elements be 
rectangular and aligned with the global X-axis, is required 
for proper stress calculation. 

4. Damping coefficients need not be included if no damping is 
desired in the dynamic analysis. 

5. The laminate input. Section III, should be symmetric as the 
displacement solution ignores bending/extensional coupling. 
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TABLE C-1. CLIP ROUTINES 


CLIP 


TIMER 

IMPIN 

LAMIN 


INPUTJ 


ELT6 

TPLATE 

STRETR, 

QTSHEL, 

QOCOS , 

TDCOS , 

TRFPRD, 

SLST, 

SLCCT 

CALBAN 


Main Program. Supervises the data input. Sets 
internal storage parameters. Supervises element 
and global stiffness formulation. 

Computes and prints elapsed CPU Time. 

Reads impact parameters. 

Reads laminate materials and orientations. Formu- 
lates laminate stiffness matrices. Computes ef- 
fective plane stress matrices for shell element 
routines . 

Reads modal data. Determines equation numbers. 
Sets impacted mode numbers to degree of freedom 
numbers. 

Sets storage for shell element routines. 

Supervises shell element formulation. 

Shell element routines. 


Finds global stiffness band width and saves elemental 
matrices for global stiffness formulation. 
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ERROR 

INL 

ADDSTF 

ADSTF2 

STEP 

STATIC 

ADDMAS 

MASSIN 

SOLSTP 

TRIFAC 

REDVK 

PR 

STRREC 

PLYSTR 

FAIL 

BSTIF 

INVRTS 
MAT PRD 
STOPPR 


Compares required storage and available storage 
Reads and processes static nodal loads. 

Forms global stiffness r mass and force matrices. 
Computes impact velocity for momentum conservation. 
Modifies existing global stiffness matrix for 
damage incorporation. 

Supervises time integration analysis. Sets in- 
ternal storage for integration analysis. 

Supervises solution for static displacement vector. 
Converts blocked mass and force vectors (from 
ADDSTF) to single, unblocked vectors. 

Reads mass vector Into core after damaged element 
reformulations. Modifies global mass and stiffness 
for impact mass separation. 

Performs time integration analysis. 

Decomposes stiffness matrix. 

Solves for displacement vector. 

Prints displacements, velocities and accelerations. 
Supervises stress recovery procedure. 

Performs composite stress analysis. 

Performs failure analysis. 

Computes properties for interlaminar shear cal- 
culations in damaged elements. 

Matrix inversion routine. 

Matrix multiplication routine. 

Terminates program execution. 


Ill 
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